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Abstract 

Inhomogeneous dynamical mean-field theory is employed to calculate the vertex-corrected elec- 
tronic charge transport for multilayered devices composed of semi-infinite metallic lead layers cou- 
pled through a strongly correlated material barrier region. The barrier region can be tuned from 
a metal to a Mott insulator through adjusting the interaction strength and the particle filling. We 
use the Falicov-Kimball model to describe the barrier region because an exact expression for the 
vertex corrections is known, allowing us to determine their effect on transport. The dc conductiv- 
ity is calculated and we find the effects of the vertex corrections are relatively small, manifesting 
themselves in a small reduction in the resistance-area product. This reduction saturates in absolute 
magnitude as the barrier layer becomes thick, as expected due to the vanishing nature of the ver- 
tex corrections in bulk. The vertex corrections have a larger relative effect on the resistance-area 
product for more metallic and thinner devices. 



I. INTRODUCTION 



As the characteristic length scales of electronic devices reach further into the nanometer 
regime, quantum-mechanical effects give rise to novel properties and new devices. The 
field of strongly correlated materials makes possible the calculation of various properties of 
nanostructures. These calculations are motivated by the significant experimental activity in 
multilayered nano-heterostructures^l New physical phenomena have been seen to emerge 
in these heterostructures that are absent in the bulk materials they are composed from, like 
the appearance of two-dimensional electron gases at the interfaces between band and Mott 
insulators^ and their low temperature superconductivity^. 

In our work, we construct strongly correlated electron nanostructures, with carefully 
controlled quantum confinement effects. These theoretical nanostructure devices, are con- 
structed from semi-infinite ballistic-metal leads that are coupled together through a barrier 
region of strongly correlated material. We investigate the effect that the barrier region inho- 
mogeneity has on the longitudinal charge transport through the multilayered nanostructures. 
Unlike previous worlP, that neglected the vertex corrections due to their vanishing nature in 
the bulkP, we investigate the effects of including the vertex corrections in charge transport. 
The vertex corrections vanish in the bulk due to the odd parity of the velocity operator 
in momentum space and the even parity of the local irreducible vertex^. The arrangement 
of the multilayers breaks the homogeneity of the system meaning this parity argument no 
longer holds for longitudinal transport. Momentum is no longer a good quantum number 
in the longitudinal direction, so one cannot change the sign of the momentum and classify 
states according to their parity. Of course, for transport parallel to the planes, we have full 
translational invariance, so the vertex corrections do vanish in dynamical mean field theory. 
Recently, using different methods from our work, the importance of the vertex corrections in 
the optical conductivity calculations of the single-band Hubbard model using the dynamical 
cluster approximation has been investigated^, showing that as frequency and doping levels 
increase, the vertex corrections play an increasingly important role. Their technique does 
not require explicit calculation of the irreducible charge vertex, but instead they directly 
calculate the optical conductivity. Additional work by V. Janis and V. PokornjP' investi- 
gated general properties of the vertex corrections to the electrical conductivity of electrons 
scattered on random impurities. They found that the sign of the vertex corrections to the 
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Drude conductivity is negative, which disagrees in sign, with the results presented in this 
paper, although the system they investigate is different from ours. Our work is concerned 
with the DC limit and investigates the effect the barrier thickness has on the vertex correc- 
tions. Due to the vanishing nature of the vertex corrections in the bulk, when the barrier 
is thick enough, it becomes bulk-like, so the vertex corrections stop having an effect, hence 
the maximal impact of the vertex corrections will be localized to the interface regions. 

The theoretical framework for this work is based on dynamical mean-field theory (DMFT) 
which allows for self-consistent calculations of the properties of strongly correlated materials. 
The DMFT formalisrrP incorporates all forms of transport (except localization effects from 
weak disorder), therefore we don't have to make any a priori assumptions about the type of 
transport in the device. We do make the assumption that the self energy is local within each 
plane, but can vary from plane to plane. For the interaction, we use the Falicov-Kimball 
modePand employ the Potthoff- No lting^ technique of a mixed basis for the Green's functions 
in solving the inhomogeneous DMFT. This involves stacking two dimensional homogeneous 
x — y planes in a longitudinal ^-direction (see Fig. [I]), with hopping allowed only between 
neighboring planes in the z-direction. We allow for inhomogeneity in the ^-direction, but 
preserve the translational invariance within each plane. This implies we can use momentum 
in the x and y directions, but we must solve the problem in real space in the z-direction. 

To begin, Fourier transform the x and y coordinates to wavevectors k x and k y , respec- 
tively, but keep the z coordinate in real space. In our formalism, we represent the z co- 
ordinate, which indicates the plane numbers, with Greek letters (a, (3, 7,...). The Green's 
function then depends only on momentum through the two-dimensional band energy, so for 
each two-dimensional band energy, we solve a quasi-one- dimensional problem that is repre- 
sented tridiagonally in real space and can be solved with the quantum zipper algorithm^. 
We iterate our many-body equations to achieve a self-consistent solution. 

In the remainder of this paper, we will work through the mathematical formalism in Sec. 
II. Section II also uses a Kubo formula for the electronic change transport including vertex 
corrections in DMFT to find the dc conductivity matrix. We initially express the polar- 
izability matrix on the imaginary-time axis, then we Fourier transform to the Matsubara 
frequencies, we take the analytic continuation to the real axis, and finally take the limit of 
the frequency as it goes to zero to find the dc-conductivity matrix and the resistance. In 
Sec. Ill, we present our numerical results for various Falicov-Kimball interaction strengths 
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FIG. 1. Schematic representation of the model, showing the stacking of homogeneous 2D layers 
forming a 3D multilayer nanostructure. The figure shows a metal-barrier-metal junction with only 
five metal layers on each side of the junction, although we can vary both the barrier and lead 
thicknesses. In the work presented here, the actual number of metal lead planes is 30 on each side 
and the number of barrier planes varies from one to fifty. 

and barrier layer thicknesses. Finally, in Sec. IV we give our conclusions and remark on 
possible future directions of research. 

II. FORMALISM 

To calculate the full electronic transport along the longitudinal direction (perpendicular 
to the planes), we use a Kubo-Greenwood-based linear- response formalism! 10 * 11 !. We operate 
in the steady state for current flow, with the charge density fixed as a function of time. A 
constant charge density implies that the charge current is conserved throughout the device 
because of the continuity equation. We also assume that temperature is a constant through- 
out the device and there is no electronic charge reconstruction 12 in our system, which is 
realistic when both the leads and the barrier are at half-filling and their chemical potentials 
match. The Kubo-Greenwood formula is based on the current-current correlation function 
of the charge current operator. The charge current operator for the ath plane is the sum of 
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all the current flowing longitudinally from the ath plane to the a + 1st plane, 




(1) 



where t aa+ \ is the hopping between the ath and a + lst plane and c ai and c ai are the creation 
and annihilation operators, respectively, for a conduction electron at lattice site % on plane 
a. Note that the index % denotes the lattice site on each plane, a square lattice for this work, 
and the planes are all aligned in registry so the longitudinal hopping is straight across each 
plane from one equivalent site to the next equivalent site, i.e. the lattice sites are those of 
a simple cubic lattice. The totally current operator is then j lon9 = } l ° n9 . 

We consider a Hamiltonian that involves a hopping term for the electrons and an in- 
teraction term for the sites within the barrier region. For the interaction in our numerical 
calculations, we employ the Falicov-Kimball modeP which involves an interaction between 
spinless conduction electrons and spinless localized electrons. When the conduction electron 
hops onto a site occupied by the localized electrons it feel a Coulomb repulsion. When this 
correlation strength is large enough, it has a Mott-like metal-insulator transition. In the 
second quantization formalism, the spinless Falicov-Kimball HamiltoniarP is, 



where £ m -j is the intraplane hopping between nearest neighbor sites on plane a, U a represents 
the interaction strength on plane a, and w a i is a classical variable that equals one if there 
is a localized particle at site i on plane a and zero if there is no localized particle at site 
i on plane a (a chemical potential fi is employed to adjust the total conduction-electron 
concentration). 

It is important to note that the following derivation does not depend on the Hamiltonian 
used and can be thought of as a general result, until we explicitly state otherwise. We add 
the perturbation, 



to the Hamiltonian, to introduce the effect of the electric field. A(i) is the vector po- 
tential, with the electric field written in the gauge where the scalar potential vanishes, 






(3) 



a 
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E a = —d t A a (t) and units chosen so c = 1 (we also choose h = 1). The electric field, vector 
potential, and charge current all are vectors in the longitudinal (z) direction only. 

To calculate the charge transport, we must start by evaluating correlation functions in 
imaginary time. To do so, we need to determine the Wick rotation for H'(t) — > H'(t). 
Since the vector potential couples to the charge current, we take the vector potential on the 
imaginary time axis to be periodic in j3 — 1/T, so it can be expressed as a Fourier series 
in the bonsonic Matsubara frequencies A a (r) = T £^ exp[iisir]A a i with iu t = 2mTl. The 
current-current correlation function for the multilayer nanostructure is defined to be 

IWM = [ dre^ (T T jt^(r)j^(0)} , (4) 

with a and (3 being planar indices (not to be confused with the (3 in the integration 
limit which is the inverse temperature) and T T denotes time ordering of operators. The 
notation (X) denotes the trace, Tr exp(— (5H)X divided by the partition function Z= 
Tr exp(—/3H), and the operators are expressed in the Heisenberg representation X(t) = 
exp(riJ) X exp(— tH), all with respect to the equilibrium Hamiltonian H . 

We calculate the current-current correlation function by taking the functional derivative 
of the appropriate expectation value with respect to the vector potential. We will use Green's 
functions to express our results. Because we need to know the dependence of the Green's 
function on A, we need to define a two-time Green's function in the presence of H + H'{t). 
The Green's function, in real space, is defined by 

G aPij (r, r'; A) — — (r T c m (r; A)c^{r'- A)) , (5) 

for imaginary time r and t' . The extra notation "; A" is used to remind us that the 
operators q q and ctg evolve according to the full Hamiltonian H + H'{t) via 0(r; A) = 
W(t, 0) O U (r, 0) with U (r, 0) being the evolution operator with respect to H + H'(t); i.e. 
d T U (r, 0) = [H+H'(t)]U(t, 0). To properly express the Green's functions for the Matsubara 
frequencies we use a double Fourier transformation 

G afi (iu n ,iu m ;A) =T [ dr [ dr'e^G^r'; A)e~^' . (6) 

Jo Jo 

We do this because the perturbation j • A(r) causes the Green's function to lose its time- 
translational invariance. Note that after taking all relevant derivatives and algebraic manip- 
ulations, we can express the Green's functions in zero field, where G a/ 3(iu n ,iu m ; A) oc 5 nm , 
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and we write G a /3(iu n ) because the limiting form of the Green's function in Eq. ^ depends 
only on r — r 1 in equilibrium as A(r) — > 0. 

To build our model, we need to solve for the local Green's function on each plane, 
which we do by employing the quantum "zipper" algorithm^, based on the Potthoff-Nolting 
formalisrrP. When solving for the local Green's functions, we use the symbol Z to represent 
a general variable in the complex plane. Typically Z takes the form of the fermionic Mat- 
subara frequencies, Z = iu n or the analytic continuation to the real axis with iu n — > u±i0 + . 
We start with the unperturbed equation of motion (EOM), for equilibrium, where A a = 0, 

^2 G~ f p(Z; fc 11 ) \{Z + n - ejjn) 5 ia + {t a _ la 5 ia 

8ya+l) ~ ^a{Z)S Ja j — S a p, (7) 
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where the two-dimensional band structure is e^ N = — 2t a [cos k x + cos/cj, fc" = (k x , k y , 0) 
is the transverse momentum, 5 a p is the Kronecker delta function, and H a (Z) is the local 
self-energy on plane a. Note that from this point on we will use the simplification that the 
hopping matrix elements are equal to t for all nearest neighbors (interplane and intraplane), 
ta+ia = ta-ia = taij = t and vanish otherwise. Since the EOM has a tridiagonal form 
with respect to the spatial component, it can be solved with the so-called renormalized 
perturbation expansion^. We solve the equation directly, for the (3 = a case via 

G ^ = \ 77m c r^ilT- ( 8 ) 

Z/ -f (J, L, a \Zj) t h \\ a -t Gaa ( Z;fc ||) t-l- Gaa (Z;kU) 1 

We create left and right recursion relations, 

L a _ n (Z; fcH) = Z + »- E a _ n (Z) - e fc „ + ^ - ^T^h ' ( 9 ) 



and 



t 2 



R a+n (Z; fc 11 ) =Z + fi- Y> a+n (Z) - e fe n + — (10) 

respectively. We start these relationships with the bulk values (n — > ±oo), which give us 

L.^Z; fcH) = Z + ^ ~ ^oo(Z) - e kU ± ^ - - _ _ _ — (n) 

and 

Roo(Z; fcH) = Z + /X ~ S 2 °° (Z) ~ efc " ± ^\Z + ^-T, 00 {Z)-e k ,Y-At\ (12) 

The signs in the previous two equations are chosen to yield an imaginary part less than zero 
for Z lying in the upper half plane, and vice versa for Z lying in the lower half plane. The 



self-energies (£±oo) vanish for the ballistic metal leads used here. Employing these recursive 
relationships, we can solve the diagonal, a — (3, and off-diagonal terms, (a ^ /?), using the 
relations 

G.„-„(g;fc") = - G 7-" + ;f^'" < 13 > 

and 

G Qa+n (Z;fcll) = -%±^|^, (14) 

which are defined for n > 0. In addition, one needs to note that the following identity holds: 
G a /3(Z;k^) = Gp a (Z;k^), by interchanging planar subscripts in the EOM, solving for the 
new recursion relationships, and showing they are the same as the above. 

We need to express the current-current correlation function in a form where we can use 
the local Green's functions for each plane. We begin by using Green's functions to evaluate 
the expectation value of the current-current correlation function in Eq. (|4]), 

SGpp+i (iu) m ,iu) n ;kW;A) 5Gp +1 p (iu m ,iu n ; k^;A) 
where 6 denotes a functional derivative. After employing the identity 



ran fell 



(15) 



G a p{iuj n ,iu m ;A) = E E G aj (iu n , iu m >] A)Gj(iu m ,, iu n r,A)G S p{iu) n >,iu m ] A), (16) 

m'n' 7<5 



we have 



n Q/3 M= E EE 



iaet 



SGj (iu m <,iu n >] k^;A) 



SA a -i 

mnm'n' fel yS 

-Gp +ll (iu m ; fc") G S/ 3 (iu n ; A; 11 )] 5 mm i5 nn >. (17) 

Note that we replaced Green's functions that have no derivative with respect to the vector 
potential by their zero field values as is done in the Kubo-Greenwood formula. 

Calculating the functional derivative of the inverse of the Green's function is easier then 
calculating the functional derivative of the Green's function itself. We find the inverse of the 
Green's function by examining the EOM for the Green's function in a field, a generalization 
of Eq. ^ to include two-time dependence and the vector potential, yielding 

y^y^^Ks^; A) [(iu m + fi - e° N ) 5 ai 5 mm > + (t5 7a _i + t8 10l+1 ) 5 mm > 

7 m' 

} = S a5 , (18) 
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with G~ 1 (icu m , iu) m i\ fell; A) being the terms inside the brackets. Taking the functional deriva- 
tive of G~ l (iu m , iu) m i\ fell; A) with respect to A a ~i, only three terms have nonzero derivative; 
resulting in 

SG~g (iu) m i, iiO n i\ fell; A) ^ 5S 7 (zu; m /, zc; n /; A) o~G 77 (iu) m i>, ioj n »\ A) 

-MaeT (M a7+ i5 a7 - t5 57 -i5 Q+ i 7 ) 8 m > +l>n >. (19) 

Evaluating the Green's functions in zero field [G (iui n i, ico n ; fc") oc <5 nm ] and noting that 
any Green's function that does not have a derivative acting on it can be replaced by it's 
(diagonal) zero-field value, the current-current correlation function becomes 

~R a p{ivi) = aVt 2 T^^ [Gpa (za; m ;fc") G a+lf}+1 (iu m+ i; fe") 

m fell 

+G/3 +la+1 (iu m ; fe") G a p (iu m+ i; fe") - G^a+i (iw m ; fe") G a/3+ i (iw m+i ; fc") 
-G/3 +la (iw m ; fc") G a+ i^ (^w; fc 11 )] 
+ iaet ^ ^ [G / 3 7 (zw m ; fc") G 7/3+ i (iu n ; fc") 

mnm"n" fell 7<5 

-Gp+n (vj m ; fc 11 ) G 7/3 («w n ; fc 11 ) r ,. ; tt " (.20) 

/J oG 77 (za; m », za; n //; ^4) M a _i 

The second term on the right hand side is the vertex correction which disappears in the bulk 
due to parity arguments, but here is nonzero. The functional derivative of S with respect 
to G is proportional to the irreducible charge vertex. 

It is worth nothing again that up until this point we have not made use of the Falicov- 
Kimball model. Hence, the formalism is independent of the model used (periodic Anderson, 
Falicov-Kimball, Hubbard, etc.). The Falicov-Kimball model is now explicitly used because 
the result for the irreducible charge vertex F is known to be^, 
• • _ 1 ££ 7 ( 

1 S 7 (iw m ) — S 7 (iw n ) 

1 7 (lW m , lW n ) - — — ; -r - r — — r<) m m"<W'> l/U-J 

1 0G 7 [lU) m », t(jJ n "', A) 1 G^lOJm) — G^ilUJn) 

in the limit of vanishing field. The delta functions on the frequencies are specific to the 
Falicov-Kimball model, other models will have more complicated formulas than what we 
develop below. 



Substituting this expression and Eq. (19) into Eq. (20) we are left with the full equation 



for the current-current correlation function on the imaginary frequency axis: 
TL a p{ivi) = e 2 a 2 t 2 T^2 \ ^2 \- G P a ( iuJm] G a+W+i ( iuj m+f, fc") 

m I fell 



fell) 

-Gp +la (*w m ; fc 11 ) G a+ip (iu m+ i; fc 11 )] 
+ ^ I ^2 ( iuJm; ^Vh ( ia W! _ G /3+i7 ( iuJ m ; fc") G 7( g (zw m+i ; fc 11 )] 

7<5 V fell 

£ 7 (iu; m ) - S 7 (zu; m+i ) 



G 7 (iu; m ) - G 7 (ia; m+ i) 



I- T\\\ ^s( tUJ m) - Esfam+l) n (■ f 

d 7 5 - > &ys(iuj m ; k 11 )——. r — — -Gs^u n ] fc 1 

4^ G s {uj m ) - G s [ioj m+ i) 



^2 \G Sa (iu m \ k ll )G a+ls (iuj m+ i; fc") - G Sa+1 (iu m ; k ll )G aS (iuj m+ i; fc") > , (22) 




x 



x 

where fc",fc",fc" are all transverse momenta in two-dimensions and the inverse denotes a 
matrix inversion with respect to the planar indices 7, 5 of the matrix defined in the brackets. 

To calculate the electronic charge transport we need the real-axis response. This is 
found by performing an analytic continuation from the imaginary to the real frequency 
axis^. We rewrite the Matsubara summations using contour integrations in the standard 
way that encircles all the simple poles corresponding to the Matsubara frequencies. The 
contour integrals have contributions at the poles of the Fermi-Dirac distribution [f(oj) = 
1/(1 + exp(/3oj))] which lie at the fermionic Matsubara frequencies, iu n = inT(2n + 1) 
(the residue of the pole is —T). The contours are then deformed to lines parallel to the 
real axis and the Green's functions are evaluated with either retarded or advanced functions 
depending on the argument and the regions of analyticity (since the functions are all analytic, 
there are no additional poles). We then set the Fermi-Dirac function at u — %v\ equal to 
f(oj) before analytical continuing the iv\ frequency to the real axis. Finally, we continue iv\ 
to v + i8 and shift the u + v — > u in relevant integrals to obtain the real axis current-current 
correlation function. 

To find the dc conductivity, we need to take the limit as the frequency v goes to zero of 
the polarization, 

a a s(0) = limRe^i^. (23) 

■u-5-0 V 

The final result is obtained after a lengthy algebraic manipulation focused on rearranging 
the a and (3 labels and explicitly writing out the real and imaginary parts, Re(Z) = ZH 2 Z * 
and Im(Z) = ^§f-, and expanding the paired terms. We also use the fact that G a p = Gp a 
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to arrive at our final expression for the de-conductivity, 



o-a/3(0) 



e 2 a 2 t 2 



du 



df(u 



J |y~] [imG^o;; feUjlmGa+^+i^; fc 1 



77 ./ \ du 
+ImG )3+ i a+ i(w; fe ll )ImG J a( g(a;; fc") - ImG^ a+ i(a;; fc^ImG^+^u;; fc 1 
-ImG^ + i a (w; feUjlmGa+i^w; fc 11 )] 
+ ]Tlm [G^(w;fcH)G 7/J+1 (w;fcll)] 

1 



fci 



x 



Im£ 7 (u;) 
Im(x 7 (u;) 



5js — 



a;;jfe") 



ImG^u;) 



Im [c^w; k ll )G a+lS (u; k 



(24) 



In numerical calculations, since the k dependence is always through e k , we replace all 
momentum summations by integrals over the 2-d DOS. Note that one can also arrive at 



Eq. (24) by using the general analytic continuation formulas in Ref. [161 



Given the conductivity matrix, we can extract the resistance of the multilayer nanos- 
tructure. This is done by relating the electric field and the expectation value of the planar 
current density in linear response and utilizing Ohm's law, giving the resistance-area per 
unit cell product of the nanostructure 



R n a 2 = y^jcr 



-1 



(25) 



Note that it is the resistance R n that is extracted from these calculations, not the resistivity. 
Since the system is inhomogeneous, we do not have the appropriate geometrical factors to 
convert resistance into the resistivity as one can do in a homogeneous system. 

The full algorithm for our work using dynamical mean-field theory can now be summa- 
rized. We start with the self-consistent iterative algorithm outlined in Ref. [4] to calculate 
the self-energies and diagonal Green's function on all planes (G aa ). We use a similar parallel 
approach exploiting the fact that the frequencies are completely independent of each other, 
to calculate the Green's functions at different frequencies on different processors. After cal- 



culating the diagonal Green's functions, we use Eqs. (13) and (14) to calculate all of the 



off-diagonal Green's functions. Once all of the Green's functions have been calculated we 



can use them in Eq. (24), noting that at half filling, the temperature dependence only enters 
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into the Fermi factor derivative, allowing for us to calculate the temperature dependence 
of cr Q( g(0) without having to recalculate the Green's functions, though all numerical results 
presented here will be at one temperature. 

III. RESULTS 

We perform our calculations at half-filling (fj, = 0, (cjcj) = 1/2, and w± = (wi) = 1/2). 
This has a number of advantages. First, because the chemical potential is the same for 
the metallic leads and the barrier, the filling remains homogeneous throughout the system 
and there is no electronic charge reconstruction. Second, the chemical potential is fixed as 
a function of temperature (/i = 0), so we don't have to recalculate the chemical potential 
as the temperature is changed. We carry out our calculations on a simple cubic lattice 
allowing for nearest neighbor hopping only (interplane and intraplane hopping are equal). 
This allows us to reduce the number of parameters that we vary for our calculations and 
focus on the physical properties. In our calculations, we also include 30 self-consistent planes 
in the metallic leads both to right and left of the barrier, which we vary in thickness between 
1 and 50 planes. 

We vary the Falicov-Kimball interaction strength (U) through its metal-insulator transi- 
tion, which for our simple cubic lattice occurs at U = 2y6 ~ 4.92. We plot the resistance- 
area product versus thickness for four different U values (Fig. [2|: U = 1 and U = 2, in 
the diffusive metal regime; U — 4, a strongly scattering, anomalous metal; U — 6, a Mott 
insulator with a small correlation-induced gap. These plots are both for vertex-corrected 
and non vertex-corrected expressions for the conductivity matrices. The plots in Fig. [2] 
show the expected behavior; for U = 1, U = 2, and {7 = 4, all show an Ohm's law lin- 
ear scaling of resistance with thickness, characteristic of diffusive metals, with a nonzero 
intercept denoting the non-vanishing contact resistance associated with the metallic leads. 
Additionally, the resistance-area product for the Mott insulator, {7 = 6 grows exponentially 
with thickness due to its tunneling behavior. Note that the vertex corrections hardly have 
any effect on the behavior in either regime. 

There are a number of numerical details that need to be discussed. We need to ensure 
that at low energies the properties of the nanostructures are accurately determined, which 
becomes increasingly difficult as we examine thicker Mott insulators and therefore requires 
more computational power. Due to these limitations, we limit our thickness for various 
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FIG. 2. (Color online) Resistance-area product for nanostructures with U = 1,2,4, and 6, and 
various thicknesses. Panel (a) is a semilogarithmic plot, while panel (b) is a linear plot. The 
temperature is T = 0.01 in both panels. Note how the correlated insulator (U = 6) has an 
exponential growth with thickness as expected^ for a tunneling process at low enough temperature, 
while the metallic cases (U = 1,2, and 4) have nearly perfect linear scaling. The values with 
(symbols) and without (lines) the vertex correction are indistinguishable on the scale of the graphs. 



interaction strengths, U = 1 up to 50 layers, U = 2 up to 30 layers, U = 4 up to 20 layers, 
and U = 6 up to 15 layers. To accurately calculate the sharp frequency dependence of 
the self-energy, we use a step size as fine as 0.001 for the frequency grid and up to half a 
million points for the integration over the two-dimensional DOS. The combination of the 
exponential growth of the resistance-area product when (7 = 6 and the inherent numerical 
accuracy we achieve makes it difficult to extract the difference between vertex corrected and 
non vertex corrected results for thicker barriers. 

In Fig. [3j we plot the difference between the vertex-corrected and non vertex-corrected 
terms to examine the effects of vertex corrections on various thicknesses. As expected the 
vertex corrections are relatively small and as the thickness increases the vertex corrections 
saturate in magnitude. This saturation is because as the thickness increases the barrier 
approaches the bulk limit, and in the bulk the vertex corrections vanish. Therefore the 
effect of vertex corrections lies primarily in the interface region and they disappear as you 
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FIG. 3. (Color online) Difference between vertex-corrected and non vertex-corrected resistance val- 
ues as a function of barrier thickness. Panel (a) shows the difference between values of the resistance 
for various interaction strengths, U, showing that the difference grows with an increase in barrier 
width until it plateaus at a saturation value that decreases as the interaction strengths increases 
towards the metal-insulator transition value. Panel (b) shows the percentage difference between 
vertex-corrected and non vertex-corrected resistance values [(R no vertex — Rvertex) / Rno vertex x 100]. 
Note how the lower the interaction strength is, the thicker the barrier must be before the saturation 
is reached, indicated by the flat region in panel (a). 

move away from interface and into the "bulk", part of the barrier, which dominates for 
thicker barriers. 



IV. CONCLUSIONS 



In this work, we used inhomogeneous DMFT to calculate the electronic dc charge trans- 
port, including the effect of vertex corrections in the Falicov-Kimball model. We used a 
Kubo- Greenwood formula to derive the dc conductivity matrix, using the exact expression 
for the vertex corrections, known for the Falicov-Kimball model. Although we only applied 



the vertex corrections to the Falicov-Kimball model, the results up to Eq. (20) hold for any 
general model. Therefore, the theory can be extended to other models if the corresponding 
irreducible vertex can be found. We can easily modify the model by adding mean-field- 
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like interactions such as Zeeman splitting for magnetic systems, and use whatever impurity 
solver is desired to calculate the local Green's functions on each plane. 

In addition to changing the general model, there are a number of transport effects and 
model modifications that could yield interesting results beyond what is presented in this 
paper. We showed that vertex corrections play a small role in dc conductivity but did 
not investigate transport outside the dc limit. We also neglected charge reconstruction^^ 
at the interface and by performing the calculations off of half-filling one could introduce 
temperature dependence to the chemical potential. Finally, by breaking the particle-hole 
symmetry one could investigate the effect of vertex-corrections in thermal transport. But in 
all cases, we expect the vertex corrections to be small for the same reasons they are small 
here. 

By varying the barrier thickness, we see the effects the vertex corrections play as we go 
from the single barrier layer limit to the bulk limit. In Figure |3j we see that the vertex 
corrections make a small relative change to the dc conductivity saturates as the barrier 
thickness increases towards the bulk. We also see that vertex corrections become relatively 
smaller as the Falicov-Kimball interaction strength increases. Hence, the effect of the vertex 
corrections are small enough that they probably do not need to be included in an inhomo- 
geneous DMFT calculation for longitudinal transport. Since the thermal transport arises 
from the Jonson-Mahan Theorem^ vertex corrections also should not affect the thermal 
transport much. 
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